% input_imex is a matlab function who gets from imex .dat file:

%  i, j, k; original grid size
%  DI ,DJ ,DK; direction block dimensions
%  PHI; An array of geostatistically distributed values is used
%for porosity
%  KI, KJ, KK; An array of geostatistically distributed values is
%used for permeability;

% by jose augusto f. ramos (zeaugustofr@gmail.com)                        %
%                                                                         %
% last updade date: 22/08/2011 

% (by now this work fine with the Tenth SPE Comparative Project .dat)

function [i,j,k,DI,DJ,DK,PHI,KI,KJ,KK,wells] = input_imex(arquivo)

fid = fopen(arquivo,'r'); %fid is file id
 
line = fgetl(fid); % read the first line of the file

%** Reservoir Description Section
%*GRID *CART
while ~any(strfind(line,'*GRID *CART')); %read the file until ...
    line = fgetl(fid);
end

mesh = line(13:28);
mesh = str2num(mesh); % ... gets the mesh size

%pre-locate new matrices 
KI = zeros(mesh(2),mesh(1),mesh(3));
% KJ = KI;
% KK = KJ; %seems that dont need this one

% line (j) - column (i) - deep (k)
[j,i,k] = size(KI);

ncells = i*j*k; % total number of cells (original size)

%** I direction block dimensions
while ~any(strfind(line,'*DI *CON'));
    line = fgetl(fid);
end

DI = line(10:14);
DI = str2num(DI);

%** J direction block dimensions
while ~any(strfind(line,'*DJ *CON'));
    line = fgetl(fid);
end

DJ = line(10:14);
DJ = str2num(DJ);

%** K direction block dimensions (layer thicknesses)
while ~any(strfind(line,'*DK *CON'));
    line = fgetl(fid);
end

DK = line(10:14);
DK = str2num(DK);

%*POR
while ~any(strfind(line,'*POR'));
    line = fgetl(fid);
end 

PHI = dlmread('input\spe_phi.inc');
PHI = reshape(PHI',[1 ncells]);
PHI = reshape(PHI,[i j k]);
 % we do this reshape and permute things to gets the permeability and
 % porosity matrix in the shape of the reservoir for smooth the upscaling
PHI = permute(PHI,[2 1 3]);

for x = 1:ncells
       if PHI(x) == 0
        PHI(x) = 0.0001;
        % "Were required some changes in the matrix of porosity, since the
        % simulator does not complete execution of this "data set" when
        % a cell with null porosity is found in the path of a well.
        
        % To solve this problem, we substitute the value of zero porosity
        %   by "0.0001", which is a value close to null." (Filho, D. S. M.)
       end
end

 
%*PERMI *ALL            ** An array of geostatistically distributed
%                       ** values is used for permeability in I direction.
while ~any(strfind(line,'*PERMI *ALL'));
    line = fgetl(fid);
end
    
KI = dlmread('input\ki.inc');
KI = reshape(KI',[1 ncells]);
KI = reshape(KI,[i j k]);

KI = permute(KI,[2 1 3]);

%*PERMJ *ALL            ** An array of geostatistically distributed
%                       ** values is used for permeability in J direction.
while ~any(strfind(line,'*PERMJ *ALL'));
    line = fgetl(fid);
end

KJ = dlmread('input\kj.inc');
KJ = reshape(KJ',[1 ncells]);
KJ = reshape(KJ,[i j k]);

KJ = permute(KJ,[2 1 3]);

%*PERMK *ALL            ** An array of geostatistically distributed
%                       ** values is used for permeability in K direction.
while ~any(strfind(line,'*PERMK *ALL'));
    line = fgetl(fid);
end

KK = dlmread('input\kk.inc');
KK = reshape(KK',[1 ncells]);
KK = reshape(KK,[i j k]);

KK = permute(KK,[2 1 3]);

% %*GROUP 'ALL-WELLS' *ATTACHTO 'FIELD'
while ~any(strfind(line,'*WELL '));
    line = fgetl(fid);
end

count = 1;
while any(strfind(line,'*WELL '));    
    %wells(:,1) = index number
    %wells(:,2) = injector = 1 / producer = 0
    %wells(:,3) = x coordinate 
    %wells(:,4) = y coordinate
    wells(count,1) = count;
    if any(strfind(line,'Injector'))
        wells(count,2) = 1;
    else
        wells(count,2) = 0;
    end
    wells(count,3:4) = str2num(line(30:36));
    count = count+1;
    line = fgetl(fid);
end
% wells = 0;

%Closes the input fid
fclose(fid);
close all;